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We present a detailed description of the recently proposed numerical renormalization group method 
for models of quantum impurities coupled to a bosonic bath. Specifically, the method is applied 
to the spin-boson model, both in the Ohmic and sub-Ohmic cases. We present various results for 
static as well as dynamic quantities and discuss details of the numerical implementation, e.g., the 
discretization of a bosonic bath with arbitrary continuous spectral density, the suitable choice of a 
finite basis in the bosonic Hilbert space, and questions of convergence w.r.t. truncation parameters. 
The method is shown to provide high-accuracy data over the whole range of model parameters 
and temperatures, which are in agreement with exact results and other numerical data from the 
literature. 
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I. INTRODUCTION 

The Numerical Renormalization Group (NRG) is 
known as a powerful tool for the investigation of quantum 
impurity problems, where a quantum system with a fi- 
nite number of internal degrees of freedom (the impurity) 
couples to an infinite system of non-interacting fermions 
with a continuous density of states (the bath) 1-5 . The 
NRG combines numerically exact diagonalization with 
the idea of the renormalization group, where progres- 
sively smaller energy scales are treated in the course of 
the calculation. NRG calculations are non-perturbative 
and - thanks to the logarithmic energy discretization - 
are able to access arbitrarily small energies and temper- 
atures. Besides providing thermodynamic quantities like 
susceptibility, entropy, and magnetization, the NRG can 
be used to compute dynamic observables directly on the 
real frequency axis. 

While the NRG was originally developed by Wilson 1 
for the Kondo model, it was later applied to a va- 
riety of more complex impurity models with one or 
more fcrmionic baths, being able to handle, e.g., two- 
channel and multi-impurity physics 6 ' 7 . As a recent ex- 
tension, impurity models with a fcrmionic bath and 
a single bosonic mode have been treated, with the 
so-called Anderson-Holstein impurity model being the 
paradigmatic example 8 . Interesting applications of the 
NRG include its use within dynamical mean-field theory 
(DMFT) 9 ' 10 . There, the electronic self-energy of a lattice 
model of correlated electrons is approximated by a local 
function in space, and the lattice model is mapped onto a 
single-impurity model supplemented by a self-consistency 
condition. Using DMFT-NRG, the Mott transition of the 
Hubbard model has been investigated in detail, both at 
zero and finite temperatures 11 ' 12 . 

The objective of this paper is an important extension 
of the NRG method, namely the application to quantum 



impurities coupled to a bosonic bath with a continuous 
spectral density (in contrast to the single boson mode in 
Ref. 8). We have recently given a short account on this 
development 13 ; the purpose here is a detailed description 
of this novel NRG application. To be specific, most of our 
presentation will focus on the spin-boson model, with the 
Hamiltonian 

H = -y cr x + ( -<r z + ^2 + yE Al ( ffl! + a i ) • 

i i 

(1) 

This model naturally arises in the description of quan- 
tum dissipative systems 14,15 : The dynamics of the two- 
state system, represented by the Pauli matrices cr XjZ , is 
governed by the competition between the tunneling term 
A and the friction term Ai(dj + at). The a, constitute a 
bath of harmonic oscillators responsible for the damping, 
characterized by the bath spectral function 

J(w) = ttJ2\^S(lu -ui) . (2) 

i 

Clearly, most interesting are gapless spectra, J(w) > 
for < uj < oj c , with ui c being a cutoff energy. In the 
infrared limit, the energy dependence of J(u>) for uj — > 
determines the system's behavior, where power-laws are 
of particular importance. Discarding high-energy details 
of the spectrum, the standard parametrization is 

J(uj) — 2irau 1 c r s uj s , < uj < oj c , s > -1 . (3) 

The case s = 1 is known as Ohmic dissipation 14 , where 
the spin-boson model has a delocalized and a local- 
ized zero-temperature phase, separated by a Kosterlitz- 
Thouless transition (for the unbiased case of e = 0). 
In the delocalized phase, realized at small dissipation 
strength a, the ground state is non-degenerate and rep- 
resents a (damped) tunneling particle. For large a, the 
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dissipation leads to a localization of the particle in one of 
the two a z eigenstates, thus the ground state is doubly 
degenerate. 

Bath spectra with exponents s > 1 (s < 1) are called 
super-Ohmic (sub-Ohmic): In the super-Ohmic case, the 
system is always delocalized with weak damping; the sub- 
Ohmic case is more involved and will be discussed be- 
low. Besides simple power-law spectra, more complicated 
bath properties can arise in a number of situations, e.g., 
structured baths, consisting of an Ohmic part and modes 
sharply peaked at certain energies, have been considered 
recently 16,17 . 

The spin-boson model has found applications in a wide 
variety of physical situations 14 ' 15 : mechanical friction, 
damping in electric circuits, decoherence of quantum os- 
cillations in qubits 18 20 , impurity moments coupled to 
bulk magnetic fluctuations 21 , and electron transfer in bi- 
ological molecules 22 ' 23 . 

Considering this wealth of applications, numerical 
methods to reliably deal with the spin-boson and re- 
lated models for all temperatures are highly desirable. 
In the past, quantum Monte-Carlo simulations 24 have 
been used, which, however, cannot work at arbitrarily 
low temperature, and cannot easily extract dynamical 
information on the real frequency axis. Density-matrix 
renormalization techniques, as employed in Ref. 25, cir- 
cumvent this problem, but are not able to resolve very 
small energy scales. In Ref. 26, the NRG with a fermionic 
bath has been used, exploiting the well-established map- 
ping of the Ohmic spin-boson model to an anisotropic 
Kondo model. Such a mapping is only valid for frequen- 
cies w<w c (which nevertheless encompasses most of the 
interesting physics), and, more seriously, is restricted to 
the Ohmic case 27 . 

In a recent paper, we have presented a formulation of 
the NRG directly for a bosonic bath, and applied it to 
the spin-boson model 13 . While we could accurately re- 
produce known results for the Ohmic case, we also found 
that the sub-Ohmic model displays two phases as well (in 
agreement with Refs. 28, 29) which are separated by a 
non-trivial quantum phase transition. Remarkably, this 
phase transition was not systematically investigated be- 
fore. We have studied the properties of the corresponding 
quantum critical points - in the phase diagram (see Fig. 1 
of Ref. 13) those form a line, parametrized by the bath 
exponent s, which terminates in the Kostcrlitz-Thouless 
transition at s = 1. Near s — 1 we could make contact 
with analytical renormalization group results, originally 
formulated by Kostcrlitz in the context of an Ising model 
with long-range l/r 1+s interaction 30 . 

The purpose of this paper is (i) to present in detail the 
implementation of the bosonic NRG method for the spin- 
boson model, (ii) to discuss various strategies to set up 
the iteration scheme for the bosonic NRG, (hi) to demon- 
strate its feasibility by studying, in particular, the case 
of Ohmic damping; and compare our data with a variety 
of results from the literature, and (iv) to discuss possible 
future applications of the bosonic NRG. The physics of 



the sub-Ohmic spin-boson model is very rich due to the 
presence of a line of boundary quantum critical points; 
a full account of the universal critical behavior, studied 
using analytical and numerical methods, will be given in 
a forthcoming publication 31 . 

The remainder of the paper is organized as follows: 
In Sec. II we introduce the formulation of the NRG 
for bosonic systems and highlight important differences 
which occur compared to the fermionic NRG. In par- 
ticular, the choice of appropriate bosonic basis states, 
which are required to accurately describe certain strong- 
coupling fixed points, is discussed, with details given in 
Appendix A and B. Section III analyzes the NRG flow 
and the low-energy fixed points, the phase boundaries 
and issues of numerical convergence as function of the 
discretization parameters. In Sec. IV we turn to thermo- 
dynamic observables calculated using the bosonic NRG, 
such as entropy and specific heat, together with their 
scaling behavior. Finally, Sec. V is devoted to dynam- 
ical quantities, where we focus on the symmetrized im- 
purity spin autocorrelation function. We close with a 
summary and discussion of applications and extensions 
of the bosonic NRG. 



II. THE BOSONIC NRG 

The bosonic NRG can be applied to a wide range 
of quantum impurity problems involving a bosonic bath 
with a continuous (and in particular gapless) spectrum. 
The following discussion of the technical details of this 
method concentrates on the spin-boson model, the first 
application of the bosonic NRG 13 . The general concepts 
are valid for the study of other bosonic impurity models 
as well. 

The purpose of this section is twofold: we want to 
discuss in detail the technical steps of the calculations in 
Ref. 13. In addition, we introduce an alternative strategy 
to set up the NRG procedure, the so-called stor-NRG, in 
contrast to the chain-NRG used in Ref. 13. The use of 
the star-NRG is related to the choice of an optimized set 
of basis states. As will be discussed in Sees. II F and 
HID, the star-NRG allows for an efficient construction 
of the NRG basis which solves the problem of the boson 
number divergence occurring in the localized phase for 
sub-Ohmic damping. 

Let us start with a form of the spin-boson model which 
is most convenient for the NRG procedure: 

l l 
H = H loc + J dsg(e)4a E + ^ J de h(e)(a E + a\) (4) 
o o 

with H\ oc = —Aa x /2 + ea z /2. In this model, g{e) 
characterizes the dispersion of a bosonic bath in a one- 
dimensional representation, with upper cutoff 1 for e. 
The coupling between the spin and the bosonic bath is 



2 



, J (CO) 




00/ CO, 



FIG. 1. Logarithmic discretization of the bath spectral 
function in intervals [A~ (n+1) , A~ n ] (n = 0,1,2,...); typical 
values of the NRG discretization parameter A as used in the 
bosonic NRG are A = 1.5 . . . 3.0. 
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given by h(e). These two energy-dependent functions are 
related to the spectral function J(w) via 



-J(x) 



de{x) 
dx 



h 2 [e(x)} (x£[0,uj c }), 



(5) 



where e(x) is the inverse function of g(x), g[s(x)] = x. 
For a given J(x), eq. (5) does not determine both g and 
h independently. Therefore, as shown below in eq. (10), 
a specific choice of h is used to simplify the calculations. 



A. Logarithmic Discretization 

The NRG procedure starts by dividing the interval 
[0,1] into intervals [A~(" +1) , A~"] (n = 0,1,2,..., see 
Fig. 1). The width of each interval is 



cL 



A" n (l - A" 1 ) 



(6) 



Within each interval we introduce a complete set of or- 
thonormal functions: 



Vv( £ ) 



1 e iu> n ps for A -(n+l) < £ < A -n 

VdZ ; (7) 

outside this interval 



(p = 0, ±1,±2, . . . and oj n = 2-n:/d n ). The operators a-p 
appearing in the Hamiltonian (4) can be represented in 
this basis: 



a e = ^2 a np ip np (e) 

np 



(8) 
(9) 
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We then choose the function h(e) to be a constant h n in 
each intervall of the logarithmic discretization: 



h(e) = h n 



1 
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FIG. 2. (a) Structure of the spin-boson model correspond- 
ing to eq. (14) as used for the chain-NRG and (b) to eq. (11) 
as used for the star-NRG; the boxes indicate the iterative 
diagonalization scheme for both cases. 

for e e [A~ (rl+1) , A~"]. With this choice, the impurity 
(the spin-operator a z ) couples to the p = component of 
the bosonic operators a np and ajj p only (the same strat- 
egy has been used in the case of a fermionic bath with 
non-constant density of states, see Ref. 32). 

The next step is to write the Hamiltonian (4) in the ba- 
sis a np and a^ p ; the p^0 components of these operators 
are still present through their coupling to the p = com- 
ponents in the free bath term. The main approximation 
of the bosonic NRG at this point is to drop this coupling, 
in close analogy to the fermionic case (see Refs. 1, 2). 
This approximation becomes exact in the limit A — > 1. 
Nevertheless, a careful check of its validity is necessary 
and will be discussed in Sec. IIIB. 

With the p^0 components completely decoupled from 
the impurity, we drop the p = index in the operators 
= o and a^ p=0 and arrive at a Hamiltonian of the form: 



^np 

with 



oo oo 

-ffioc + V £„at a n + —±= V in (a„ + a\) , (11) 

n=0 2 V T n=Q 



tn=ln 2 dxJ( X )x, ri= dxJ(x). (12) 

JA-(»+i)u c JA-("+i)w c 

The label H s is introduced to distinguish this "star"- 
Hamiltonian from the "chain" -Hamiltonian H c (see 
eq. (14) below). The and the 7„ can be easily eval- 
uated for a bath spectral function of the form given in 
eq. (3): 

, a + 1 1 - A-( s + 2 ) A 
U = ~ ~ : — ^-rrW^A 



s + 2 1 - A-( s +!) 
27ra 



7„ = 



L w 2 ^ _ A -(-+l)^ A -n(-+l) 



(13) 



The structure of this Hamiltonian is sketched in 
Fig. 2b: the impurity spin couples linearly to all the 
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bosonic degrees of freedom a n , in a very similar way as 
in the original Hamiltonian (1). The bath spectral func- 
tion for model (11) is discrete, consisting of 5-peaks at 
energies £„ with weight oc 7^. Each bosonic degree of 
freedom of this star-Hamiltonian is a representative of 
the continuous spectrum of bosonic degrees of freedom 
in the intervals [A~ (n+1) , A~™]. 



B. Chain-NRG vs. Star-NRG 

Starting from the model (11), there are two possi- 
ble ways to set up a numerical renormalization group 
procedure. The first one (which we call "chain" -NRG 
in the following) uses the transformation of the star- 
Hamiltonian (11) to a semi-infinite chain: 



H c — H\ oc - 

00 

+ [ e nb f n b n + t n (blb n+1 + 6^ +1 6 n ) 



n=0 



(14) 



with 770 = J dx J(x) . The spin now couples to the first 
site of the bosonic chain only (see Fig. 2a) and the re- 
maining part of the chain is characterized by on-site en- 
ergies e„ and hopping parameters t n , in analogy to the 
fcrmionic NRG. The parameters e„ and t n can be calcu- 
lated numerically from a given spectral function J(w), as 
discussed in detail in Appendix A. 

Such a mapping from a star-Hamiltonian on a semi- 
infinite chain form is exact. It has been used in all ap- 
plications of the fermionic NRG since the original work 
of Wilson 1 . Its generalization to the bosonic NRG is 
straightforward and has been employed in Ref. 13 (see 
also Ref. 33). 

The structure of the Hamiltonian (14) is sketched in 
Fig. 2a. The boxes indicate the NRG strategy used in this 
case: in the first step, a cluster containing the impurity 
plus the first bath site is diagonalized. In each subse- 
quent step, the cluster is enlarged by one additional site 
and the new cluster is diagonalized using the information 
obtained in the previous step. 

The second possibility (which we call "star" -NRG in 
the following) is to use the Hamiltonian (11) directly 
for the iterative diagonalization. The general idea is 
sketched in Fig. 2b: again, the first step of the renor- 
malization group procedure involves the diagonalization 
of a cluster containing the impurity plus the first bath 
site. The following renormalization group steps, how- 
ever, are completely different to the chain-NRG as each 
new bosonic site does not couple to the previously added 
site but to the impurity instead. 

The suggestion to use such a star-NRG for the investi- 
gation of bosonic impurity models, such as the spin-boson 
model, raises a couple of questions: 



(2) Is the star-NRG of any advantage as compared 
to the chain-NRG (apart from the simplification 
that we do not have to calculate the e„ and t n of 
cq. (14))? 

(3) Why has such a star-NRG not been used in the 
fermionic case? 

The answers to questions No. (1) and (2) will be given 
further below. Let us first discuss question No. (3) in 
more detail. A fermionic star-NRG for, say, the Kondo 
model would start from a Hamiltonian similar to eq. (11). 
The important difference in the fcrmionic case is that 
the logarithmic discretization has to be performed for 
both positive and negative frequencies. As a conse- 
quence, there are two sets of bath operators in the star- 
Hamiltonian, one for positive and one for negative fre- 
quencies: 

oc OO 



<r,n— 



CT,n— 



For a hybridization function close to particle-hole sym- 
metry we have w £~ . This means that at each renor- 
malization group step one has to add two fermionic sites 
(the alternative to add f^} n+ first and then f^\_, or vice 
versa, suffers from violating particle-hole symmetry, if 
present). The Hilbert space therefore increases by a fac- 
tor of 16 in each step. It is much more convenient to 
first map the star-Hamiltonian to a chain form similar to 
eq. (14). In this form, only one site has to be added in 
each renormalization group step. 

Whether such a fermionic star-NRG is of any advan- 
tage is not clear. It might be useful for extreme asym- 
metric cases, but for the cases which are usually of inter- 
est the chain-NRG already works very well and is much 
easier to implement. 

Coming back to the bosonic NRG, there does not seem 
to be an a priori preference for either star- or chain-NRG 
because the structure of the bosonic bath is extremely 
asymmetric from the outset (restricted to positive fre- 
quencies only). To address the possible advantages of the 
star-NRG, we first have to give more details of how the 
bosonic NRG is implemented (for both star- and chain- 
NRG). 



C. Iterative Diagonalization and Choice of Bosonic 
Basis States 

The star-Hamiltonian H = H s (11) and the chain- 
Hamiltonian H = H c (14) can be written as a series 
of Hamiltonians Hm (N > 0) equal to H in the limit 

N -> 00: 



H = lim A~ N H N . 

N^oo 



(16) 



(1) Does the star-NRG work at all? 



The Hn for the star-Hamiltonian are given by 
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Hn.s — A w 



N 



N 



+ 77= ln ( a " + a " 

n=0 V n=0 

(17) 

and for the chain-Hamiltonian by 



H 



N,c 



A N 

N 



N-l 

Y 6 nbtK + Y f ™ ( fo n 5 «+l + b l+l b r 



71 = 



n=0 



(18) 



In this notation, both Hq_ s and i?o,c correspond to a two- 
site Hamiltonian with only the first site of the star or 
chain coupled to the spin. 

Two successive Hamiltonians are related by the follow- 
ing renormalization group transformations: 



H 



N+l,s 



AH 



N,s 



+ A N+1 



Cjv+iajv+iOw+i + 2^^ N+1 ( aN + 1 + a lv+i) 



(19) 



and 



H 



N+l,c 



AH 



N.i 



■ A N+1 



(20) 



The factor A N in eqs. (17) and (18) enables the direct 
comparison of the low-frequency spectra of subsequent 
Hamiltonians and, in particular, the discussion of fixed 
points as in Sec. III. In contrast to the fermionic case, 
the factor is A N instead of A N / 2 because the energies £„ 
in the star-Hamiltonian and the e n and t n in the chain- 
Hamiltonian are falling off as A - ™, instead of the t n oc 
A - ™/ 2 in the fermionic case. (This implies that a bosonic 
NRG calculation with discretization parameter A and a 
particle-hole symmetric fermionic one with A 2 will have 
comparable energy resolution.) Note that, in the sub- 
Ohmic spin-boson case, the 7„ are falling off slower than 
A~ n . Nevertheless, the factor A~™ is the appropriate one 
for the low-energy spectra as shown in Sec. HE. 

The sequences of Hamiltonians (17) and (18) are solved 
by iterative diagonalization. In the first step, the H are 
diagonalized in a basis formed by the product states of 
tr z -eigenstates \a) and a suitable basis for the first bath 
site (we will describe below what we mean with 'suitable 
basis'). We have to introduce a cutoff Abo already for 
this basis, but this is usually not a serious restriction as 
we can use fairly large values of A^o ~ 500 (in contrast to 
the much lower values of A^ for the following iterations) . 

Given the eigenstates |r)jv 01 



H N \r) N = E N {r)\r) N 



1,...AT S 



(21) 



with N 8 the dimension of Hn, we can construct a basis 
of H N+1 : 



\r;s) 



N+l 



\r) N ®\s{N + 1)) 



(22) 



with \s(N + 1)) a suitable basis for the added site. In 
setting up the basis \s(N + 1)) we are faced with two 
problems not present in the fermionic case: 

(1) The numerical approach restricts the number of ba- 
sis states one can take into account to a maximum 
number N\> ~ 10 — 14. The validity of this approx- 
imation has to be checked carefully. 

(2) A criterion for a suitable selection of ATb basis states 
out of the infinitely many states of the added site 
has to be found. 

A general criterion for an "optimal" basis (for a given 
A^b) can be formulated as following: find a set of bo- 
son states \s(N + 1)) which give the best description of 
the lowest-lying many-particle states of iJjv+i (see also 
Ref. 34). In a variational sense, this corresponds to find- 
ing states ^(A - -!- 1)) which give the lowest many-particle 
energies for a whole set of energy levels (see also Fig. 4 
below) . This is certainly not a rigorous statement and we 
have not yet developed a general algorithm to set up such 
an optimal basis. Instead we select one of the two sets of 
basis states optimized for the two stable fixed points of 
the spin-boson model: the A^ eigenstates of 6jv +1 6jv+i 
(or (i i N+1 aN + i) with lowest eigenvalues as an optimal ba- 
sis for the delocalized fixed point (Sec. II D) and displaced 
oscillators as optimal basis for the localized fixed point 
(Sec. HE). 

Before continuing let us point out that are no symme- 
tries in the Hamiltonians Hn. s and Hn^ c (at least for the 
interesting case of finite a and A). This is in contrast to 
the fermionic case 1,2 , where we can use, for example, the 
total spin and particle number as quantum numbers to 
significantly reduce the size of the Hamiltonian matrices 
(which would be of size (4A f s ) 2 in the absence of symme- 
tries). Consequently, in the bosonic NRG for the spin- 
boson model there is only one matrix of size (A^A^) 2 to 
be diagonalized in each renormalization group step. This 
results in a much simpler structure of the NRG program, 
but limits the values of N s to 100-200. 



D. Optimal Basis for the Delocalized Fixed Point 

Let us start from the a = limit of the spin-boson 
model in which two-level system and bosonic degrees of 
freedom are completely decoupled; for finite A, the spin 
oscillations are undamped and the system is in the delo- 
calized phase from the outset. 

The Hamiltonian in the original formulation (1) then 
takes the form: 



A t 
H = -—a x +2_^Wia\L 



(23) 
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For simplicity, the bias e is set to zero. The star- 
Hamiltonian in the a = limit has the same structure: 

A x , 

H s = -jcr x + 2_^ ■ (24) 

n 

From this structure, it is clear that the eigenstates 
of a^ N+1 a N+ i with lowest eigenvalues form the optimal 
basis: 

\s(N + l)) = {\n N+1 )} , (25) 

with 

a N+i a N+i\n-N+i) = n\n N+ i) , n = 0,1, . . . N h - 1 . 

(26) 

The reason is simply that here the many-particle energies 
are given by the sum of the single-particle energies £ n . 

The situation is similar in the chain-NRG, where the 
a = limit reads: 

A °° 

H c = - —a x + ^2 [tnbib-a + t n {b ] n b n+ i + • 

n=0 

(27) 

Here we choose for the basis \s(N + 1)) the states 
{\n N+ i}} with 

b N+i b N+i\n N +i) = n\n N+ i) , n = 0,1, . . . N h - 1 . 

(28) 

The difference to the basis for the star-Hamiltonian is 
that the |njv+i) are not eigenstates of the full bosonic 
part in eq. (27). But in contrast to the case of a > 0, the 
Hamiltonian (27) conserves the total number of bosons; 
the many-particle states with the lowest energies are then 
given by those states which are constructed from the 
single-particle states with the smallest boson numbers, 
independent of whether a diagonal basis is chosen or not. 

In our previous implementation of the bosonic NRG 13 
we used the basis (28). This is a suitable choice only if 
the many particle states of J/jv+i with lowest energies are 
indeed constructed from states with small boson number 
- in other words, if the average values of the boson num- 
bers (&jy +1 &Ar + i) are small. This is the case when the 
system is close to the delocalized and the quantum crit- 
ical fixed points. However, the boson number diverges 
when the system flows to the localized fixed point for 
s < 1 as discussed below. 



E. Optimal Basis for the Localized Fixed Point 
(Displaced Oscillators) 

Here we consider the spin-boson model with zero tun- 
neling amplitude, A = 0. In this case, oscillations be- 
tween | |) and | |) are absent and the system is in the 
localized phase from the outset. 



The Hamiltonian in the original formulation (1) then 
takes the form: 

i i 

For simplicity, the bias e is set to zero. As the bath de- 
grees of freedom now couple to a static spin, the Hamil- 
tonian can be decomposed in two sectors ffj for a z = +1 
and for a z = —1: 

H 1 =Y ff *t ' H ^ = w i a \ a i + y ( a i + 4) > ( 30 ) 

i 

(Hi accordingly). In each sector, we now have indepen- 
dent bosonic degrees of freedom which can be written 
as: 

Hit = Uiajai , (31) 
(dropping a constant term) with 

a, =ai + 6i , 6 t = A . (32) 

The quantities 9i can be viewed as an effective (dimen- 
sionless) coupling between impurity and bath mode i. 
Apparently, this transformation corresponds to a dis- 
placement of the oscillators ai by the value +9i for the 
T-sector and — 6i for the .[-sector. The displacements do 
not change the energies Wj. This means that the whole 
many-particle spectrum of the bosonic bath is identical to 
the one for the uncoupled bath except for the additional 
two-fold degeneracy corresponding to the two sectors f 
and J.. 

Note that for the original spin-boson model (1) the 
LUi and Ai are not specified independently, only the bath 
spectral function J(u>) is given; therefore we cannot give 
explicit expressions for the 0i for eq. (1). 

The star-Hamiltonian eq. (11) for A = (and e = 0) 
takes a form similar to eq. (29). Again we have two 
sectors with 

x 1 ■ 

H s \ = + -z~j= In {a n + ai) , (33) 

n=0 n=0 

(H s i accordingly). Using the same reasoning as above, 
we can now write 

oc 

Hsi = Y ' ( 34 ) 

with 

a n = a n + 9 n , 6 n = . (35) 

2V 7r ?n 

The values of j n and are given in eq. (13) so we obtain 

9 n oc A^ 1 -*)/ 2 oc C ( r 1)/2 . (36) 
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Written in terms of energy oj we therefore have: 

9(lo) oc u^- 1 ^ 2 . (37) 

This result is rather interesting: for sub-Ohmic baths, 
s < 1, the shift 9 n grows exponentially with n. However, 
in the super-Ohmic case the shift goes to zero in the low- 
energy limit (n — > oo), and it is energy-independent for 
the Ohmic case. Technically, the coupling to the impurity 
can be viewed as a relevant (irrelevant) perturbation of 
the discretized spin-boson model for s < 1 (s > 1) and as 
a marginal perturbation in the Ohmic case. Thus, in the 
Ohmic and super-Ohmic case the effective coupling #(w) 
does not diverge as u> — > even in the extreme localized 
case of A = 0. Therefore, numerical problems associated 
with a diverging effective coupling are only expected in 
the sub-Ohmic case. 

Coming back to the iterative numerical diagonalization 
of the star-Hamiltonian, it is now clear that a simple basis 
as in (25) can be far from the optimal choice. If we stay in 
the original basis constructed from the lowest eigenstates 
of a^cin, we need more and more basis states to describe 
the lowest eigenstates of the displaced oscillators. 

On the other hand, it is clear how to construct the 
optimal basis for \s(N +1)) at least for the A = case. 
For the sectors f / J. we simply take oscillator states with 
displacements +0n+i/ — 9n+i- As we need a single basis 
for both sectors, these states have to be orthogonalized 
first; this will be discussed in more detail in Appendix B. 

The displaced oscillator states can also be used to di- 
agonalize the chain-Hamiltonian (14) for A = 0. For 
a given iteration number N, the H Nfi for the f-scctor 
reads: 



H 



N.cl 



N-l 

+ e « fe « 6 « + + b i+i b n) 



Introducing displaced oscillators 

K =b n + 9 n (N) 
we again have a diagonal form 

N 



ffjv,cT = A N e n blb n . 



• (38) 



(39) 



(40) 



n=0 



The displacements 9 n (N) can be calculated numerically 
for any given set of {e„} and {t n }. For fixed N they 
show the same qualitative behavior as the n for the star- 
Hamiltonian: 



\9 n (N)\ oc A nl - 1 ^ s ^ 2 . 



(41) 



It turns out, however, that the 9 n (N) depend on both 
n and N with significant deviations from the exponen- 
tial form for n close to N . This has important con- 
sequences for the use of displaced oscillators as basis 



states in the chain-NRG. Let us assume that we used 
±9n(N) to construct the basis for -ffjv.c. Adding the site 
N + 1 introduces a significant change in the displacement 
6n(N) — > 9jy(N + 1). One possible solution to this prob- 
lem is to anticipate the coupling to the site N + 1 by 
adding a static displacement term, which is subtracted 
again in the next step. Such an approach gives correct 
results for the chain-NRG when we set A — 0. We did 
not, however, succeed in implementing the displaced os- 
cillator idea for the general case of finite A in the chain- 
NRG. So far, this strategy only works for the star-NRG 
as described in the following subsection. 



F. General Strategy of the Bosonic NRG 

In the preceding subsections we have described various 
options of how to set up the bosonic NRG. We have intro- 
duced both a star- and a chain-representation of the spin- 
boson model and we discussed two possibilities for choos- 
ing a basis for the added site: eigenstates of ^jv+i^jv+i 

(or aj v+1 OAr+i) as in Sec. IID or displaced oscillators as 
in Sec. HE. 

Now we want to discuss how we actually proceed with 
the bosonic NRG: how do we decide between the different 
options described above? 

As a starting point we choose the chain-NRG using 
eigenstates of b N+1 b^ + i (the basis denoted by |njv+i)) as 
the simplest possible basis. This approach has been used 
for all the results shown in Ref. 13. From the discussion 
in Sees. II D and II E we anticipate that this choice of the 
basis is reasonable for 

• all parameters in the super-Ohmic case, 

• the Ohmic case provided the coupling a is not too 
large, 

• and the sub-Ohmic case provided the system is 
close to the delocalized fixed point. 

On the other hand, it is clear that there will be prob- 
lems when the system flows to the localized fixed point 
in the sub-Ohmic case. The situation in the crossover 
regions and close to the quantum critical points needs 
to be checked numerically: it turns out that the critical 
fixed points for all < s < 1 can be reached using the 
\n>N+i) basis. 

There is a simple criterion to decide when the ba- 
sis |nj\r+i) is sufficient. Consider the expectation value 
njv+i(-/Vb) = (fojv + i&Ar+i) for the cluster after adding the 
site N + 1, calculated for a temperature of the order of 
the level spacing at this NRG step. This quantity can 
be obtained numerically up to values of N\> ~ 14. If the 
lowest eigenvalues of b N+1 b^ + i are a good choice for de- 
scribing the lowest eigenstates of H^+i, then nN + i(Nb) 
should be small and rapidly saturate with increasing A^. 
If, on the other hand, we identify that njv+i(-/Vb) does not 
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2 4 6 8 10 12 14 



N b 

FIG. 3. Dependence of the expectation value 
njv + i = (&jy +1 &jv+i) on the number of basis states Nb for 
a chain-NRG calculation. The values of ri]v+i(A r b) quickly 
saturate for a — 0.05 < a c and a c = 0.06113 whereas no 
saturation is observed for a — 0.1 > a c . (Parameters are 
s = 0.6, N = 20, A = 2.0, N s = 60, A = 0.01). 

saturate but increases with Nb, then we certainly have to 
abandon the basis |rijv+i) and use a different 'optimized' 
basis. 

This behavior is shown in Fig. 3 where we show results 
from the chain-NRG for a sub-Ohmic bath (s = 0.6), 
A = 0.01 and three values of a in the vicinity of the 
quantum phase transition. For a < a c and a = a c we 
indeed find a rapid saturation of HN+i(Nb) whereas no 
saturation (at least up to Nb = 14) is observed for a > 
a c . 

The behavior of riN +1 (Nb) for a > a c can be easily 
understood from the discussion of Sec. II E: as the system 
is flowing to the localized fixed point corresponding to 
the effective A approaching zero, we have to use properly 
displaced oscillators as a basis. The increase of tin+i (Nb) 
just means that we need more and more states in the 
undisplaced basis to describe the lowest eigenstates of 
Hn+i- 

In this case, the use of displaced oscillators as intro- 
duced in Sec. HE is much more appropriate. Note, how- 
ever, that the shifts 9 n cq. (35) can only be defined from 
the outset for the A = case. For any finite A, the sys- 
tem evolves according to the iterative diagonalization. If 
the system turns out to flow to the localized fixed point, 
we have to use effective displacements 6 n to set up the 
basis. These displacements have to be extracted numeri- 
cally from the renormalization group calculation and are 
different (for finite A) from the 9 n given in cq. (35). 

Figure 4 describes the general strategy to determine 
the optimal values of the displacements. The low-energy 
spectrum of ifjv+i is calculated for a whole set of 8- 
values. According to the discussion in Sec. II C, we 
identify the optimal 6 as the one which gives the low- 
est eigenenergies in iJjv+i- This value is indicated by the 
vertical line in Fig. 4, which shows results for the sub- 
Ohmic case and parameters close to the localized fixed 
point. There is a plateau in the energy levels close to the 
optimal value which means that a slight variation of the 
9 affects the lowest energies only very weakly. Note that 



uf -30 




e 

FIG. 4. Dependence of the energies of the lowest eigen- 
states of Hn+i on the displacement 9 used for constructing 
the basis for the new degree of freedom added in each itera- 
tion step. All levels shown here have their minimum at the 
same value 9 = 9* (indicated by the vertical line) which is 
the optimal value for setting up the basis. The parameters 
for this calculation are: s = 0.2, A = 1.0, a = 0.25 > a c . 

E n (6) = E n (—8), therefore a maximum at 6 = 0. The 
corresponding figure for parameters close to the delocal- 
ized fixed point (not shown here) gives a minimum of the 
many-particle levels at 9 = 0. For further details of this 
procedure, see Appendix B. 

The data of Fig. 4 are calculated using the star-NRG 
formulation. Although a similar figure can be generated 
using the chain-NRG, we are facing the (so far unsolved) 
problem discussed in HE: adding a new site changes the 
optimal displacements for the previous iterations. For 
this reason, all the results in this paper using a basis 
of displaced oscillators are calculated within a star-NRG 
representation. 

G. Diagonalization and Truncation 

To conclude Sec. II, let us briefly discuss the remain- 
ing technical steps necessary to complete the iterative 
diagonalization. For a given basis, we first set up the 
Hamiltonian matrices 

H N+1 (rs, r's') = N+i(r; s\H N+1 \r'\ s') N +i ■ (42) 

For both chain- and star-formulation of the NRG, the 
matrices can be written as a sum of three parts: 

H N+1 (rs, r's') = H$ +1 + H$ +1 + H^ +1 , (43) 

with 

#w+i( rs > r ' s ') = A N+i(r;s\H N \r';s') N+1 

= AE N (r)S rr ,S ss , , (44) 

for both chain- and star-formulation and 
#$-i,s(™,rV) 

= A N+1 £ N+1 N+i(r; s|aj vr+1 aj\r+i|r / ; s') N+1 

= A N+1 Z N+1 6 rr , (s(N+l)\J N+1 a N+1 \s'(N + l)) , (45) 
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(with operators a replaced by b for -ffjy^ c ). 

The third term takes the following form for the star- 
NRG 

ff$. liB (rs,rV) 



_ A w+i 7w+i Ar+1 ^ r; s \ az [a N+l + a\f +1 )\r'; s') N+1 
N(r\<Jz\r') N 



2V¥ 
A iv+i 7iv+i 



x(s(A + l)| aA r +1 + | S '(AT + 1)) , 
and for the chain-NRG 



(46) 



H^ +1 , c (rs,r's') 

= A N+1 t N N +i(r; s\b f N b N+1 + h.c.\r'\ s') N+1 
= A N+ H NN {r\tf N \r') N {s{N + l)\b N+1 \ S '(N + 1)) 
+h.c. . (47) 

All matrix elements of the form (s(A + 1)| . . . \s'(N + 1)) 
can be further simplified once the basis \s(N + 1)) is 
given. Similar to the fermionic case, the matrix element 
N(r\b' N \r')N appearing in the chain-NRG eq. (47) can 
be written in terms of the unitary matrices necessary 
to diagonalize H jy. The matrix elements Ar(r|(7 z |r')jv i n 
eq. (46), however, have to be calculated iteratively. (The 
technical details are very similar to the fermionic case, 
see Refs. 1, 2). 

With N s the dimension of Hm and Ab the number 
of basis states in \s(N + 1)), we then arrive at a single 
(N s ■ A b ) x (N s ■ A b ) matrix for H N+1 (rs, r's'). This ma- 
trix can be diagonalized using standard routines. From 
this we obtain the unitary matrices UN+i(rs,f) and the 
spectrum of eigenenergies i?jv+i(r) so that 

H N+1 \f) N+1 = E N+1 (f)\f) N+1 , f = 1, . . .N s ■ N h . 

(48) 

In contrast to the fermionic case, no symmetries can be 
taken into account to separate the matrix Hpf + i(rs, r's') 
into smaller sub-matrices. 

The dimension of -ffjv+i now has to be reduced from 
N s ■ Ab to N s to allow for a numerical calculation with 
computation time growing only linearly with N. This is 
achieved with the usual truncation scheme where only the 
lowest N s eigenstates of Hm+i are kept (for the fermionic 
case see Refs. 1, 2). These states form the basis |r)jv+i 
for the next step and the iteration continues. 

The calculation of correlation functions, such as the 
spin-spin correlation function C(u) in Sec. V, requires 
the calculation of additional matrix elements N(r\A\r')iy. 
For more details see Sec. V. 



III. FLOW AND FIXED POINTS 

The iterative numerical diagonalization of the spin- 
boson model as described in the previous section gives 
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FIG. 5. Flow diagrams calculated with the chain-NRG for 
the parameters s = 1 (Ohmic case), iv c = 1, e = 0, A = 0.01, 
and a = 0.6 in (a) and a = 1.4 in (b). The NRG parameters 
are iV s = 100, iVb = 8, and A = 2.0. 



a sequence of many-particle levels E^{r) (r = 1, . . . N s ). 
Due to the logarithmic discretization, these energies fall 
off as E]y(r) oc A~ N . NRG flow diagrams can then be 
constructed by plotting A n En(t) versus iteration num- 
ber TV. 

In this section we focus on those issues which can be 
directly inferred from the NRG flow diagrams: the ap- 
pearance of fixed points, the crossover between different 
fixed points at finite energy or temperature, and quantum 
phase transitions between the fixed points. Subsections 
III A-III C deal with the Ohmic spin-boson model; here 
we also address the issue of convergence. In subsection 
III D we investigate those features connected to the flow 
of energy levels which are specific for the sub-Ohmic case. 

All results are calculated for cutoff energy oj c = 1 and 
bias e — 0; we employ NRG parameter values of A = 
1.8 - 3.2, A b0 = 100, A b = 4 - 14, A s = 30 - 120. 



A. Fixed Points 

Let us first concentrate on results from the chain-NRG 
for the Ohmic case, s = 1, and various values of A and 
a. 

Fig. 5 shows two NRG flow diagrams for A = 0.01 
and two values for the coupling: a = 0.6 in Fig. 5a 
and a = 1.4 in Fig. 5b. In these diagrams, the rescaled 
many-particle energies A N E]y(r) are plotted versus the 
iteration number A with the ground state energy sub- 
tracted. Another difference to the fermionic case (apart 
from the different prefactor A^ instead of A N / 2 ) is the 
absence of an even-odd effect: in the fermionic case, the 
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many-particle spectrum usually oscillates between two 
sets of energy levels (so that it is more appropriate to 
speak of a limit cycle than of a fixed point). Plotting the 
many-particle spectrum either for even or for odd itera- 
tion numbers only then gives the flow diagrams as shown 
in, for example, Refs. 1, 2. 

In our bosonic NRG calculations, we can follow the 
flow typically up to N = 60 (corresponding to T w 10~ 20 
for A — 2.0), then we observe an unphysical runaway 
which is due to the accumulation of numerical errors in 
the course of the iteration. As the runaway scale depends 
on the numerical precision used in the program code, it 
can be shifted to lower temperatures if needed. 

The flow diagrams of Fig. 5 show the existence of two 
different fixed points: the delocalized fixed point for small 
a (see Fig. 5a, N > 20) and the localized fixed point for 
large a (see Fig. 5b, N > 6). These two fixed points are 
stable and the quantum phase transition between them 
is discussed further below. 

If the value of A is small enough (as in Fig. 5a) the sys- 
tem is close to the localized fixed point in an intermedi- 
ate range (4 < A^ < 8 in Fig. 5a) even for a values below 
the critical coupling a c . This has direct consequences for 
thermodynamic properties in the corresponding tempera- 
ture range (see, for example, Fig. 14). However, the vicin- 
ity to the localized fixed point does not imply localization 
in the sense that a system initially prepared with the im- 
purity spin in one specified direction remains in this spin 
state under time evolution. For any finite temperature, 
thermal excitations destroy localization (see Ref. 14). 

In Fig. 5a, we also observe a crossover from the lo- 
calized to the delocalized fixed point which takes place 
at A" w 10 — 20. The corresponding crossover scale, T* 
(which - in the Ohmic case - is equivalent to the renor- 
malizcd tunnel splitting A r up to a prefactor), will be 
discussed in subsection IIIC. 

The flow diagram of Fig. 5a is similar to the one ob- 
tained in Ref. 26 (Fig. 1 in Ref. 26), where the mapping 
of the spin-boson model to the anisotropic Kondo model 
was employed. The structure of the many-particle levels, 
however, cannot be directly compared as they reflect the 
type of bath used in the NRG approach (bosonic in our 
case, fermionic in Ref. 26). 

The spectrum of the delocalized fixed point (Fig. 5a 
for A" > 20) is identical to the spectrum of a spin-boson 
model with zero coupling between spin and bosons (a = 
0). The Hjy for the chain-NRG then take the form: 



10.0 



Hn,c — A 



N 



Hi oc + ^2 e nblb ri 



n=0 



JV-1 



+ 



'n+l + K +l 



71 = 



b n ) 



(49) 



In this Hamiltonian, impurity and bath degrees of free- 
dom are completely decoupled and can be diagonalized 
separately. The spectrum of the impurity part (ifi oc = 
—Aa x /2) is non-degenerate. The bath part is that of a 
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FIG. 6. Comparison between the fixed-point spectra for 
the delocalized fixed point calculated with the chain-NRG 
(circles) for various values of Nh and the fixed point spectra 
constructed from the single particle levels co n in eq. (50) (solid 
lines) for a selection of states En(v). The parameters are 
s = 1, A = 1.0, and a < a c . The NRG parameters are 
N s = fOO, and A = 2.0. 



free chain of bosons with AT + 1 sites which can be diag- 
onalized exactly: 

N JV-1 N 

E e nb\\b n + E t n (bibn+i + &t + i&n) = E u n b f n b n . 



n=0 



n=0 



n=0 



(50) 



Figure 6 shows a comparison between the fixed point 
spectra for the delocalized fixed point calculated with 
the chain-NRG (circles) and the fixed-point spectra con- 
structed from the single particle levels w n in eq. (50) 
(solid lines). The NRG data are calculated for differ- 
ent ATb- The agreement is very good for the first few 
excitations already for Ab ~ 6, while a larger value of 
A^b is required to correctly reproduce the excitations at 
higher energies. 

While the delocalized fixed point is reached for a 
smaller than a critical a c (A), the system is in the lo- 
calized phase for all a > a c (A). The localized phase 
is characterized by a (renormalized) tunneling amplitude 
A r = and a two-fold degenerate ground state. In the 
language of the (perturbative) renormalization group 14,35 
the localized phase corresponds to a line of fixed points, 
parametrized by a. Interestingly, the fixed-point value of 
a does not influence the eigenenergies of the many-body 
fixed-point Hamiltonian, but only its eigenstates, see the 
discussion in Sec. II. Thus the NRG level spectrum in 
the entire localized phase is identical to the one for the 
delocalized fixed point, apart from an additional two- fold 
degeneracy of all many-particle levels. This feature can 
be clearly seen in Fig. 5. [Of course, the approach to the 
localized fixed point depends on the particular value of 
a, consequently the NRG flow on intermediate scales will 
be different for different a > a c (A).] 
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B. Critical Coupling and Convergence 

The results shown in Fig. 5 indicate the well-known 
transition between the localized and delocalized fixed 
points at a critical a c (A) 14,15 . Due to the Kosterlitz- 
Thouless nature of this transition, the fixed point at 
a = a c (A) is not a new fixed point, but belongs to the 
localized phase instead. 

On approaching the transition from the delocalized 
side, we find, as expected, that the crossover scale van- 
ishes as 38 InT* oc l/(a c - a), see Fig. 11 below. We 
use this dependence to determine the value of a c from 
numerically calculated data for T*(a) via a non-linear 
fit. [Note that on the localized side of the transition a 
low-energy scale only shows up in the flow towards the 
fixed point, i.e., in corrections to the fixed-point values 
of observables; thus the critical a c (A) is easier obtained 
via extrapolation from the delocalized side.] 

As already discussed in Ref. 13 (see Fig. 2 in Rcf. 13), 
the value of a c also depends on the NRG parameters A, 
Ab, and N B . Figures 7a-c show the characteristic depen- 
dence for the Ohmic case, s = 1, and two values of A 
in Fig. 7c. Keeping A fixed, we observe a rapid conver- 
gence of a c with increasing N s (Fig. 7a) and N b (Fig. 7b). 
Note that we did not observe a transition to the delocal- 
ized phase for N b < 4, even for very large values of a. 
As expected from the iterative diagonalization scheme, 
the values of N s and Ab necessary for convergence in- 
crease with decreasing A (see Fig. 7c). The converged 
data for a c (A) show a linear A dependence in the range 
1.8 < A < 3, with a deviation of about 15% at A = 2 
from the extrapolated A — > 1 value. 

We find that the slope in a c (A) is independent of A 
which is connected to fact that the logarithmic discretiza- 
tion systematically underestimates the spectral weight 
contained in J(u) (for a discussion of this point in the 
fermionic case, see eq. (5.42) in Ref. 2; for the soft-gap 
Anderson model see Fig. 4 in Ref. 36). 

The extrapolated values a c (A, A — > 1) for the Ohmic 
case are summarized in Fig. 8. In the limit of small A, 
the NRG result is in good agreement with the well es- 
tablished value a c (s = 1, A — > 0) = 1. We estimate the 
error in a c to be of the order of 0.02 which is due to the 
various extrapolations just described. The solid line in 
Fig. 8 shows a linear fit to the numerical data which gives 
a c (A) = 0.99 + 0.53 A. This is consistent with the RG 
result a c = 1 + 0(A/uj c ) 14 . 
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FIG. 7. Dependence of the critical coupling a c on the NRG 
parameters A B , Ab, and A for the Ohmic case; (a) dependence 
on A B for fixed Ab = 10; (b) dependence on Ab for fixed 
A s = 100 (A = 2.0 for both (a) and (b)); (c) A-dependence 
of a c for two values of A, and various NRG parameters Ab 
and A s . The dashed lines are linear fits to the Ab = 8 and 
A s = 100 data in the range 1.8 < A < 3. 
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FIG. 8. Dependence of the extrapolated value a c (A — > 1) 
on the parameter A for the Ohmic case s = 1. The crosses are 
the numerical data and the solid line is a linear fit which gives 
a c (A) = 0.99 + 0.53A. The NRG parameters are A s = 100 
and A b = 8. 



C. Scaling 

We expect to observe scaling behavior in all physical 
properties for fixed A and a — ► a c (A) and for fixed a 
and A — > 0. Such a scaling can already be identified on 
the level of the flow of the many-particle energies. An 
example is shown in Fig. 9 for fixed a = 0.6 and various 
values of A. In this way we can easily determine the 



crossover scale T* for the crossover from the localized 
to the delocalized fixed point (there is only a single low- 
energy scale): 

T* = const, x A~ N * , (51) 

where we define N* as the value of N where the first 
excited state reaches the value En = 0.3. Note that 
a change of this (arbitrary) value can be absorbed in a 
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FIG. 9. Scaling of the flow of the many-particle levels 
EN{r) for fixed a — 0.6, s = 1, and various values of A. 
The NRG parameters are N s = 100, N b = 8, and A = 2.0. 
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FIG. 10. Dependence of the crossover temperature T* on 
A for s = 1 and fixed values of a. The exponents in T* oc A x 
are x » 1.55 for a = 0.4, x « 2.15 for a = 0.6, and x w 3.49 
for a = 0.8. The NRG parameters are N s = 100, N b = 8, and 
A = 2.0. 

change of the prefactor in eq. (51); this reflects the fact 
that a temperature scale can only be defined up to a 
constant prefactor anyway. 

In the scaling regime, the dependence of T* on a and 
A is given by 14,38 



T* oc A 1/(Qc ~ a) 



(52) 



As shown in Figs. 10 and 11, the NRG results are in 
agreement with eq. (52). 



D. Flow for Sub-Ohmic Baths 

As already mentioned in Sec. II F, the chain-NRG with 
a basis of undisplaced oscillators as in eq. (28) is sufficient 
for the Ohmic and super-Ohmic case. Let us now turn to 
the sub-Ohmic case where we expect problems with the 
chain-NRG when the system is flowing to the localized 
fixed point. Figure 12 shows a typical flow diagram of the 
many-particle energies, calculated with the star-NRG for 
s = 0.8 and a couple of a values close to the quantum 
critical point a c = 0.40294. 

In contrast to the Ohmic case, the transition in the 
sub-Ohmic case is characterized by a new fixed point, the 
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FIG. 11. Dependence of the crossover temperature T* on 
a for s = 1 and fixed values of A (data for A = IO -3 and 
A = 10 -4 same as in Fig. 4b of Ref. 13). The values for the 
critical coupling are a c = 1.162 for A = IO -2 , a c = 1.150 
for A = 10~ 3 , and a c = 1.147 for A = 10~ 4 . The NRG 
parameters are = 100, At, = 8, and A = 2.0. 
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FIG. 12. (Color online) Flow diagram of the lowest lying 
many-particle energies calculated with the star-NRG for the 
sub-Ohmic case (s = 0.8, A = 0.1), using displaced oscillators 
as optimized basis. The critical value is a c = 0.40294. The 
NRG parameters are N B = 80, N b = 8, and A = 2.0. 



quantum critical fixed point, with a level structure which 
is different from both the localized and the delocalized 
fixed points. For any a ^ a c there is a finite crossover 
scale T* for the crossover to the localized fixed point (for 
a > a c ) and to the delocalized fixed point (for a < a c ). 
The crossover scale can be defined in a similar way as in 
Sec. Ill C. A further analysis of the dependence of T* on 
\a — a c \ gives the critical exponents. Their s-dependence 
has been shown already in Fig. 5a of Ref. 13. A detailed 
investigation of the critical properties of the sub-Ohmic 
spin-boson model will appear elsewhere. 

Here we focus on the level structure of the localized 
and delocalized fixed point in Fig. 12. Both fixed points 
have exactly the same level structure apart from an ad- 
ditional two-fold degeneracy of all levels of the localized 
fixed point. This is evident from Fig. 12 (levels for a > a c 
and a < a c converge to the same spectrum) and also 
follows from the discussion of Sec. HE. However, the 
proper description of the localized fixed point can only 
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N 

FIG. 13. Dependence of the expectation 
value tl at = {fcjyfejv) on the iteration number N for s = 0.8, 
A = 0.01, a = 0.2 < a c (squares), a = a c = 0.21488785 (dia- 
monds), and a = 0.4 > a c (circles), and two different values 
of N b , calculated with the star-NRG. 

be achieved using an optimized basis with displacements 
calculated as discussed in Sec. II F and Appendix B. Us- 
ing a basis of undisplaced oscillators (9 = 0) leads to an 
incorrect level structure. This can be seen in the upper 
right-hand panel of Fig. 3 in Ref. 13 (s = 0.6, a > a c ) 
where the basis (28) was used. The resulting hxcd point 
levels are therefore not the same as the one for a < a c 
in the upper left-hand panel of Fig. 3 in Ref. 13. 

As mentioned in Sec. II F, we did not yet succeed 
to implement the concept of displaced oscillators in the 
chain-NRG, so the proper description of the localized 
fixed point for s < 1 is presently only possible with the 
star-NRG. Fortunately, the problems of the chain-NRG 
only show up when the flow is approaching the localized 
fixed point. We can therefore safely extract all the criti- 
cal properties such as critical exponents from the chain- 
NRG, as has been done in Ref. 13. 

On the other hand, the use of a basis of displaced os- 
cillators within the star-NRG solves the problem of the 
boson- number divergence (see Sec. II F). This is illus- 
trated in Fig. 13 where the dependence of the expecta- 
tion value njv = (bjyfeAr) is shown for three values of a 
(a = 0.2 < a c , a = a c = 0.21488785, and a = 0.4 > a c ) 
and two values of A b . For all values of a we observe a 
rapid convergence with N^, similar to the convergence 
shown for a < a c and a — a c in Fig. 3. The difference 
here is that the data converge with A b also for a > a c 
which can not be achieved by using the basis (25), see 
Fig. 3. Furthermore, the expectation value tin diverges 
exponentially with N for a > a c , as expected from the 
discussion in Sec. HE. A diverging boson number itself 
is therefore not a problem for the bosonic NRG, provided 
a proper optimized basis is chosen. 

Finally, a few words on the limitations of the star- 
NRG. Whereas the localized fixed point is described cor- 
rectly, the star-NRG seems to fail in other respects: the 
low-energy flow to the delocalized fixed point appears 
incorrect, and critical exponents of the quantum phase 
transition deviate from the chain-NRG results (and from 



analytically known values). We do not yet fully under- 
stand this problem, but it might be connected to trun- 
cation errors which affect the star-NRG in a completely 
different way as the chain-NRG. (The idea is that the 
truncation somehow affects the character of the impurity 
operator to which the added bosonic site couples in each 
step.) The precise characterization of this problem and 
its possible solution are left for future studies. 

IV. THERMODYNAMIC QUANTITIES 

In this section, we describe how thermodynamic quan- 
tities can be extracted from the flow of many-particle lev- 
els Ejj(r) which are calculated with the bosonic NRG. 
Starting from the Epf(r) there is no difference (from a 
technical point of view) between the fermionic and the 
bosonic case (for the fermionic case see, for example, 
Refs. 2, 37). Nevertheless, for completeness we include 
a brief discussion of the technical details here. We show 
results for the impurity contribution to the entropy and 
the specific heat in the Ohmic case (using the chain-NRG 
with basis (28)). The Ohmic case has been studied in de- 
tail in Refs. 26, 39 (for earlier work on thermodynamic 
properties see Refs. 14, 40, 41). The agreement with the 
results from Refs. 26, 39 is excellent which again confirms 
the reliability of the bosonic NRG for the investigation 
of quantum impurity models involving a bosonic bath. A 
few comments on thermodynamic properties in the sub- 
Ohmic case are given at the end of this section. 

Consider the spectrum of many-particle energies Ei of 
a discretized version of the spin-boson model (not neces- 
sarily the discretized Hamiltonians (11) and (14) used in 
the bosonic NRG). The grand canonical partition func- 
tion, Z = Tre-PiH-iiN)^ reduces t0 

Z = -£e-f> E ', (53) 

i 

as the chemical potential fi is set to zero (we are inter- 
ested in gapless spectral functions J(oj)). Free energy F 
and entropy S are then given by 

F=-T\nZ and S = . (54) 

(We set fce = 1.) The impurity contribution to the en- 
tropy is 

Simp = S — So , (55) 

where S is the entropy of the full system and So the 
entropy of the system without impurity. 

Before we discuss the full temperature dependence of 
Simp(T), let us focus on the value of Si mp at the local- 
ized and delocalized fixed points: 5i mP; L and S'imp.D- It 
is well known that <Si mPi L = In 2 and Si mPl D = 26,39 , 
but it might not be obvious that these values can be di- 
rectly extracted from the many-particle spectra at the 
fixed points. 
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In Sec. Ill A we already showed that the fixed point 
spectrum of the delocalized fixed point is the same as 
the one of a free bosonic chain, which is nothing else but 
the system without impurity. This implies that for the 
delocalized fixed point 

Ei = E ifi + AE , (56) 

with Ei (Ei t o) the many-particle energies of the system 
with (without) impurity and AE a constant shift inde- 
pendent of i. It is clear that this equation cannot hold 
for all levels, it is only valid for energies sufficiently below 
the crossover scale to the fixed point. 

Equation (56) directly leads to the proof of <Si mPi D = 0: 
we have Z^> — exp[— f3AE]Z , and from this Fu = F + 
AE. The energy shift drops out in the derivative so that 
Sb = So and the impurity contribution to the entropy at 
the delocalized fixed point is given by Si mPi D = 0. 

In a similar way one can easily prove that Si mPi L = In 2: 
in this case we have 

Z L = 2 e~ PEl , Ei = E ifi + AE , (57) 

i 

with the factor of 2 due to the additional double degen- 
eracy of all many-particle levels at the localized fixed 
point. This gives Zl = 2 exp[— (3AE]Zq, and from this 
F L = -T In 2 + F + AE and S L = In 2 + S , correspond- 
ing to iSj mPj L = In 2. From this discussion it follows that 
Sim P ,L = In 2 and Si mp ,D = independent of the expo- 
nent s in the spectral function J(uS). 

For any finite A and a, the values Si mp ,L = In 2 and 
Sim P ,D = are strictly valid only in the limit T — > 0. 
Note that a proper definition of these zero-point entropies 
requires the correct order of limits: the thermodynamic 
limit has to be taken before the limit T — > 0. With the 
order of the limits reversed, the zero-point entropy would 
be equal to lnc? g , with d g the degeneracy of the ground 
state. Although this happens to give the same values for 
Sim P ,L and .Sim^D in the case studied here, this equiva- 
lence is not generally valid. (This can be seen, for ex- 
ample, in the NRG calculations for the single-impurity 
Anderson model 2 where the degeneracy of the ground 
state oscillates between 1 for even and 4 for odd itera- 
tions when the system approaches the fixed point of a 
screened spin, which has 5i mp = 0. Also, any non-trivial 
quantum critical fixed point is expected to have a residual 
entropy which is not lnrf g with integer d g .) 

The impurity contribution to the entropy is close to 
a fixed point value also when the system is close to this 
fixed point in an intermediate range of the flow diagram. 
From Fig. 5a we can therefore immediately see that the 
temperature dependence of <Si mp (T) contains a crossover 
from a high-temperature value Simp w 5i mPj L = In 2 to 
the low-temperature value Si mp (T — > 0) = Simp.D = 0; 
provided the flow is to the delocalized fixed point. The 
detailed behavior of «Si mp (T) in the crossover region re- 
quires a numerical calculation as described below. 



In the bosonic NRG, we do not have access to the full 
spectrum of many-particle energies Ei as used in eq. (53). 
Instead, the iterative procedure results in a sequence of 
many-particle energies E^{r) with iteration number N 
and r — 1, . . . N s . According to the discussion in Refs. 1, 
2, each of the sets of many-particle energies is assumed 
to be a good description of the system for a certain tem- 
perature Tjv with 

T N = xoj c A~ n , (58) 

with x a dimensionless constant of the order of 1, chosen 
such that T/v lies within the spectrum Ex{r). 

For each iteration step N, the partition function is 
calculated for the temperature T^: 

Zn = Y, e ~ EN{r)/TN ■ ( 59 ) 

r 

In addition, the internal energy at iteration step N for 
the temperature is calculated as 

E N = 1 Ly2E N (r)e- E "W T " . (60) 

r 

This is the information we have available for the numer- 
ical calculation of thermodynamic properties. 

One possibility to proceed is to calculate the free en- 
ergy Fm = —Tn In Zn for each iteration step, and from 
this the entropy S = — §fJ via a discrete differentiation. 
This procedure has been shown to give good results in 
the fermionic case (see, for example, Ref. 42). It re- 
quires, however, a precise calculation of the difference 
of the ground state energies between subsequent steps; 
this appears to introduce some errors in the calculations 
within the bosonic NRG. (In general, the bosonic NRG is 
less accurate in the calculation of thermodynamic prop- 
erties as compared to the fermionic NRG because we can- 
not keep as many states as in the fermionic case.) 

Therefore, we use an alternative approach in which the 
entropy Sn at iteration step N for the temperature Tjv 
is calculated via 

S N = ^+\nZ N . (61) 

J-N 

This approach avoids the discrete differentiation, and 
does not require the knowledge of the ground state ener- 
gies. 

Let us now discuss the results for entropy and spe- 
cific heat calculated with the bosonic NRG using the 
method just described. Figure 14 shows the temperature 
dependence of the impurity contribution to the entropy, 
<Simp(T), for a = 1/3, s = 1 (Ohmic case), and vari- 
ous values of A. We observe a crossover from the high- 
temperature value <Si mp = In 2 to the low-temperature 
value Si mp = at a crossover scale T* , which is the same 
as the one introduced in Sec. IIIC. The crossover scale 
decreases with decreasing A in agreement with eq. (52). 
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FIG. 14. Temperature dependence of the impurity contri- 
bution to the entropy, 5 l i mp (T), for a = 1/3, s = 1 (Ohmic 
case), and various values of A. 




FIG. 15. (a) Scaling curves of the impurity contribution to 
the entropy, S'imp(T), for s = 1 (Ohmic case), and various 
values of a; (b) Scaling curves of the impurity contribution to 
the specific heat, C\m P {T)/(T/T*), for the same parameters 
as in (a). 



Note the similarity of Fig. 14 to Fig. 9 for the scaling 
of the energy levels, a similarity which is simply due to 
the relation between Si mp (T) and the flow of the many- 
particle levels. 

As briefly mentioned in Sec. Ill A, the vicinity to the 
localized fixed point for early iterations (which results 
in the high-temperature value Si mp (T) w In 2) does not 
imply localization. The value of S[ mp (T) for high tem- 
peratures is due to the fact that for temperatures T> A 
both states of the two-state system contribute equally to 
the thermodynamics. Note also the similarity to Si mp (T) 
in the Kondo model: there the high-temperature phase 
is that of a local moment with both spin f and J. con- 
figurations contributing to the entropy (a temperature 
dependence of Si mp (T) as in Fig. 14 might therefore ap- 
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FIG. 16. Temperature dependence of the impurity contri- 
bution to the entropy, Si mp (T), in the sub-Ohmic case for 
various values of a and s = 0.8 (main panel) and various val- 
ues of s (inset). The coupling a is below a c so that the flow is 
to the delocalized phase for all parameters in this figure. Lines 
with symbols in the inset are data from the bosonic NRG and 
solid lines are fits assuming a power-law, Si mp (T) oc T s . 



pear more natural in the Kondo model but, of course, it 
is also valid here). 

The scaling behavior of Si mp (T) for fixed a = 1/3 
and various A is obvious and is shown in Fig. 15a to- 
gether with the scaling curves for a = 1/5, 1/4, and 1/2. 
The agreement with the exact results from the Bcthc 
Ansatz calculations in Ref. 39 is very good (see Fig. 7a 
in Ref. 39), in particular for the a-dependence of the 
scaling curves. 

The temperature dependence of the specific heat, 
C imp (T), is calculated via C imp (T)/T = dS imp (T)/dT. 
Here we cannot avoid the discrete differentiation of 
Si mp (T). The scaling of S- lmp (T) implies a scaling of 
Cim P (T)/T as shown in Fig. 15b. This figure is also very 
similar to previous calculations (see Fig. 2 in Ref. 26 from 
the NRG via mapping to the anisotropic Kondo model 
and Fig. 7b in Ref. 39 using the Bcthc Ansatz), and we 
find the same characteristic features here: a linear spe- 
cific heat C oc T for low temperatures, a peak in C/T 
at T w T* for small dissipation a < 0.3 in contrast to 
the monotonous decrease of C/T for large dissipation 
a > 0.3, and a characteristic crossing point of all the 
C/T scaling curves. 

Similar to the NRG calculations in Ref. 26, the thermo- 
dynamic quantities can only be calculated on a discrete 
mesh of temperatures given by eq. (58). This strongly 
limits the resolution of the peak in C/T for a < 0.3, in 
contrast to the Bethe Ansatz calculations of Ref. 39. 

The physics in the sub-Ohmic case is much richer, due 
to the appearance of a line of quantum critical points 13 . 
This is reflected in the behavior of the entropy and the 
specific heat. For the results of Si mp (T) and Ci mp (T) 
close to the quantum critical points we refer to a sub- 
sequent publication. Here we focus on the flow to the 
delocalized phase. 
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Figure 16 shows the temperature dependence of the im- 
purity contribution to the entropy, S , ; mp (T), in the sub- 
Ohmic case, s — 0.8, for various values of a below the 
critical value a c rts 0.125. For a close to a c we observe 
a two stage quenching of the entropy of the free moment 
(the quantum critical point has a non-trivial zero-point 
entropy of 5 qcp (r — ► 0) w 0.6 for s = 0.8). As expected, 
the temperature scale for the crossover to the delocalized 
fixed point increases with the distance from the criti- 
cal point. The low-temperature behavior of Si mp (T) for 
a < a c is given by Si mp (T) oc T s which can be seen more 
clearly in the inset of Fig. 16 where Si mp (T) is plotted for 
various values of s. This behavior is in agreement with 
the calculations of Ref. 40, where C(T) oc T s was found 
for the slightly asymmetric (e ^ 0) sub-Ohmic spin-boson 
model. [While the finite e turns the quantum phase tran- 
sition into a smooth crossover, it does not influence the 
qualitative low-energy behavior in the delocalized phase.] 

The data in Fig. 16 are calculated with the chain-NRG. 
The results from the star-NRG look similar (they give, in 
particular the correct value Si mp (T — ► 0) = In 2 if the flow 
is to the localized phase). We observe, however, a low- 
temperature behavior for S lmp {T) which is different from 
the correct form, <Si mp (T) oc T s . As briefly mentioned in 
Sec. HID, the reason for this failure of the star-NRG is 
presently not clear but probably due to truncation errors. 

Despite these deficiencies, the bosonic NRG is a reli- 
able tool for the calculation of thermodynamic properties 
in a wide range of parameters and the comparison with 
well-established results is very promising. Thermody- 
namic properties in the quantum critical region will be 
discussed in a separate publication. 



V. DYNAMIC QUANTITIES 

The calculation of dynamic properties is straightfor- 
ward within the bosonic NRG and proceeds in a very 
similar way (from a technical point of view) as in the 
fermionic case. The typical problems such as the com- 
bination of information from different iteration steps 
and the broadening of the discrete spectra have been 
discussed already in the literature (see, for example, 
Refs. 43-46, 12) and need not be repeated here. 



A. Dynamical Spin Correlations 

One important dynamic quantity of interest in the 
spin-boson model is the spin-spin correlation function 
(spin autocorrelation function) 



_L [ + °° 
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'C(t) dt 



(62) 



with C(t) = ^{[(T z (t),a z ] + ). We only consider equilib- 
rium correlation functions, in general for finite temper- 



atures, but the focus here is on T = so that the ex- 
pectation value (. . .) has to be taken with respect to the 
ground state. 

For a discrete Hamiltonian, the spin-spin correlation 
function at T = can be written as 



CH = i^|(0| ( T z |n)| 2 (5( W + eo-e„) , u > 0, 



(63) 



with C(ui) = C(—ui). Note that with the above defini- 
tion of C(t), the quantity C(w) is purely real and related 
to the imaginary part of the spin-susceptibility x( w ) y i a 
C(u) = f 7r|ImxM| (see also eq. (3.96) in Ref. 14). 

Due to the truncation in the course of the iterative di- 
agonalization, we cannot calculate C{oj) simultaneously 
for all energy scales. Instead, the correlation function 
is calculated for each cluster of length N (which gives 
information on energy scales of the order of A~ N ) and 
this information has to be added up properly. Finally, 
the discrete spectrum has to be broadened which results 
in continuous curves for C(uj) as shown, for example, in 
Fig. 17. 

These technical issues are dealt with using the ap- 
proach described in Ref. 12; to broaden the spectra, we 
use a Gaussian on a logarithmic scale (see eq. (8) in 
Ref. 12) with broadening parameter b = 0.7. 

We also define the correlation function S(w) as 



S(w) 



2C(w) 



The static spin-susceptibility \ is defined as 



X = 2- 



de 



, (a z ) = (0\a z \0) . 



It is related to C(u) via 



X = 4 
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and, using eq. (63), can be written in the form 

l<ok»| 2 



(64) 



(65) 



(66) 



(67) 



Here, we calculate the susceptibility according to eq. (67). 

The bosonic NRG allows the calculation of dynamic 
properties in a wide range of frequencies so that the func- 
tional dependence of C on the frequency ui (such as a 
power-law C(u) oc uj s ) can be easily extracted. However, 
the exponent of the calculated C(uu) has a deviation from 
the expected value s (for the flow to the delocalized fixed 
point) of about 2%. To extract the correct prefactor of 
C(uj) (which we use to compare with the exact results at 
the Toulouse point and to check the Shiba relation), we 
need to redefine the quantity S(iv) as 



2C(u) 



(68) 



16 




FIG. 17. Spin-spin correlation function C(u>) calculated 
for the Ohmic case, s = 1, A = 0.01, and various values 
of a < a c close to the transition. For small frequencies, 
C{u>) oc u>, whereas for higher frequencies we observe a di- 
vergence, C(ui) oc with logarithmic corrections. 
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FIG. 18. Spin-spin correlation function S(u>) = 2C(ui)/ui 
at the Toulouse point a — i for the Ohmic case s = 1 and 
various values of A. The inset shows the scaling of these 
curves (S(w)/S(0) plotted versus w/A r ) together with the 
exact result (thick dashed line). 



where 5 is the exponent fitted to C(uj) in the small fre- 
quency regime. 

In Fig. 17, C(lu) is shown for the Ohmic case and a set 
of a values close to the critical a c . The spin-spin corre- 
lation function shows the expected power-law behavior, 
C{lo) oc uj, in the low-frequency regime uj < T* . In the 
limit of a — ► a c , C(uj) shows a divergence for uj > T* , 
C(uj) oc uj' 1 , with logarithmic corrections. 

In Fig. 18, S(uj) is plotted at the Toulouse point 
(a = i) of the Ohmic spin-boson model for several val- 
ues of A. At this point, the Ohmic spin-boson model 
is exactly solvable, as discussed in Rcf. 15. In the in- 
set of Fig. 18, all the curves are rescaled onto one curve 
with a renormalized tunneling amplitude A r . Here, A r 
is defined as 



i r ,NRG 



Xnrg 



(69) 



where Xc is the exact susceptibility and A ri6 the exact 
renormalized tunnelling amplitude at the Toulouse point: 
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FIG. 19. NRG results for Xnrg calculated at the Toulouse 
point a = \ for the Ohmic case s = 1 and various values of A 
and A (circles: A = 0.0125, squares: A = 0.025, diamonds: 
A = 0.05, triangles: A = 0.1). Dashed lines show the exact 
values Xc = 16/ (tt 2 A 2 ) and thin solid lines are fits to the 
numerical results. 



A Pi e = 7rA 2 /2 and XeA r , c = 8/n. The quantity xnrg is 
the susceptibility calculated from the NRG, eq. (67). The 
comparison of the result from the bosonic NRG with the 
exact rescaled S(uj)/S(0) shows good agreement (see the 
inset of Fig. 18). The exact result is given by 15 : 



ln(l + 4x 2 ) 2arctan(2x) 



+ ■ 



, (70) 



S(uj) _ 1 
5(0) ~ 8(1 + a; 2 ) 

with x = 

The NRG results for xnrg deviate significantly from 
the exact value Xe — 16/(-7r 2 A 2 ). However, as shown 
in Fig. 19, this deviation is entirely due to discretization 
effects and the extrapolation A — ► 1 shows almost perfect 
agreement with the exact result. Note that the exact 
value for Xe has been obtained for a soft cut-off in the 
bath spectral function, J(uj) = 2Traujexp(—uj/uj c ). To 
allow for a comparison, the logarithmic discretization has 
to be performed for the same soft cut-off (we introduce 
a high-energy hard cut-off at uj = 15w c ). 

The scaling behavior of S(uj)/S(0) for fixed a and dif- 
ferent values of A is shown in Fig. 20. For this we need 
to identify an energy scale T* as in Sec. Ill C. There are, 
as usual, various possibilities to define the energy scale: 
the position of the peak in C(uj), uj p , the quantity 1/x, 
and the T* as defined in eq. (51). Obviously, we have 
T* oc uj p oc 1/x oc A r , and we choose A r = 8/(7rx) for 
the energy scale in Fig. 20. The scaling curves shown in 
Fig. 20 are in good agreement with the ones calculated 
in Ref. 47; in particular, we find that the coherent peak 
in S(uj) disappears when a is larger than a* w 0.3. 

In our notation, the Shiba relation reads 41 



2a (|) 2 = S(0). 



(71) 



Table I shows the results from the bosonic NRG for the 
Ohmic case and various values of a and A. The parame- 
ter S is the exponent defined in eq. (68). We find that the 
Shiba relation is fulfilled within an error of about 10%. 
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FIG. 20. Scaling spectra for the spin-spin correlation func- 
tion S{u) = 2C(uj)/lo s for various values of a in the Ohmic 
case s = 1. 

TABLE I. Results from the bosonic NRG for the 
Shiba-relation in the Ohmic case for various values of a and 
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B. Order Parameter 

In the localized phase, which corresponds to the or- 
dered phase of the 1/r 2 Ising model, it is straightforward 
to define an order parameter, to, corresponding to the 
magnetization of the Ising model. In the language of 
the spin-boson model, it corresponds to the static, i.e., 
ui = 0, part of the spin autocorrelation function; in the 
language of the anisotropic Kondo model this just mea- 
sures the prefactor of the Curie part of the local sus- 
ceptibility, i.e., the unscreened fraction of the impurity 
moment. The Kostcrlitz-Thouless nature of the transi- 
tion implies a jump of the order parameter at the phase 
transition. 

We extract this order parameter from the S(uj) contri- 
bution to C{uj). The identification of such a <5-peak in 
the spectrum of C(u>) requires some extra care. The spec- 
trum calculated with the NRG consists of (5-peaks only, 
which have to be broadened suitably to give spectra as 
shown, for example, in Fig. 17. Therefore, one has to de- 
cide whether a <5-peak in the spectrum belongs to the con- 
tinuum or whether it survives as a (5-peak in the thermo- 
dynamic limit 48 . The procedure is illustrated in Fig. 21. 
Let us first note that the matrix element |(0|<t z |0)| 2 van- 



10" 

10" 1 

10" 2 

10" 3 

10" 4 
1.2 



-— a=0.7 
(/ 83 
--- a=0.9 
-- a=0.96 
--- a=1.2 




FIG. 21. (a) Energy of the first excited state -B(l) versus 
iteration number for s = 1.0, A = 0.01, and various values of 
a. (b) matrix element |(0|ct z |1}| 2 versus iteration number for 
the same parameters as in (a). 



ishes for all parameters of Fig. 21. We therefore plot the 
matrix element |(0|<7 2 |1)| 2 in Fig. 21b together with the 
energy E(l) in Fig. 21a. We observe that for a > a c 
the energy E(l) vanishes faster than A"^ with increas- 
ing iteration number TV, whereas the matrix element ap- 
proaches a constant, |(0|ct z |1)| 2 — ► const w 1. In the 
thermodynamic limit, this gives the 5-peak at a; = 0, 
with the weight given by the matrix element |(0|er z |l)| 2 
which corresponds to the order parameter to. On the 
other hand, for a < a c the energy E(l) is proportional 
to A~ N , and the corresponding <5-peak is therefore inter- 
preted as being part of the continuum. 

These arguments result in an order parameter to (a) 
which is zero for a < a c and jumps to a finite value for 
a > a c . In the sub-Ohmic case, the order parameter 
shows power-law behavior near the quantum phase tran- 
sition, which will be discussed in detail elsewhere. As 
an aside, we note that the order parameter can also be 
extracted from the Curie part of the static local suscep- 
tibility X (T). 



VI. CONCLUSIONS 

In this paper we have discussed a generalization of Wil- 
son's NRG technique to quantum impurity problems with 
a bosonic bath. Focussing on the application to the spin- 
boson model, we have shown that this novel method pro- 
vides reliable results for both static and dynamic quanti- 
ties in the whole range of model parameters and tempera- 
tures. For the case of Ohmic damping, we have compared 
our data to existing results in the literature and found 
good agreement. The bosonic NRG is able to reproduce 
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the expected scaling behavior as function of temperature 
or frequency. For sub-Ohmic damping, there is a line of 
continuous boundary quantum phase transitions for all 
< s < 1, with exponents varying as function of s; de- 
tails of the associated quantum critical behavior will be 
discussed in a forthcoming paper (see also Refs. 13, 31). 

We have outlined several details of the numerical im- 
plementation of the bosonic NRG. Two general strategies 
were discussed, termed chain-NRG and star-NRG: both 
use a sequence of boson states with exponentially de- 
creasing energy scales, but in the chain-NRG the bath 
states form a chain and the impurity couples to the first 
chain site only, whereas in the star-NRG the impurity is 
coupled to all bath sites which are not connected to each 
other. The advantages and disadvantages of both meth- 
ods were discussed in detail, together with the important 
issue of the optimal choice of a basis set of bosonic states 
at each bath site. This problem is inherent to the bosonic 
NRG, as the infinite Hilbcrt space has to be truncated, 
and specific solutions have to be found for the problem at 
hand. We have argued that in the Ohmic, super-Ohmic 
and sub-Ohmic cases of the spin-boson model (except for 
the flow to the localized fixed point in the latter case) the 
basis formed by the lowest boson number eigenstates is 
sufficient, and all fixed points are properly captured in 
the NRG. Most of the results in this paper were obtained 
with this basis choice using the chain-NRG method; we 
have given a detailed account on convergence issues with 
respect to truncation and discretization parameters. In 
the sub-Ohmic case, the boson numbers diverge in the 
localized regime in the low-energy limit, and a different 
basis choice is needed. We have described suitable basis 
states using displaced harmonic oscillators, which solve 
the problem for the star-NRG. Open numerical issues in- 
clude a reliable implementation of the displaccd-oscillator 
basis for the chain-NRG, a more accurate calculation of 
dynamic quantities, as well as the numerical stability for 
very long iterations, i.e., very small energy scales. 

The bosonic NRG can be easily generalized to im- 
purities with multiple bosonic baths or both fermionic 
and bosonic baths. This is the subject of current 
work and will allow the study of large classes of im- 
purity models, e.g., so-called Bose Kondo 21 and Bosc- 
Fermi Kondo models 49 . These models are known to dis- 
play intermediate-coupling fixed points associated with 
universal local- moment fluctuations. The Bose-Fermi 
Kondo model arises in the context of extended dynami- 
cal mean-field theory (EDMFT) 50 , where a lattice model 
is mapped onto an impurity model with a fermionic 
bath (representing conduction electrons) and a bosonic 
bath (representing bulk spin fluctuations). The quan- 
tum phase transition appearing in the Bose-Fermi Kondo 
model has been proposed to describe local quantum crit- 
ical behavior in EDMFT, which may be relevant to the 
physics of certain heavy-fermion quantum phase transi- 
tions. However, a full numerical solution of the EDMFT 
equations at T — has not been presented to date, due 
to the lack of suitable impurity solvers. A version of the 



bosonic NRG may help to overcome this difficulty. 

Other applications of the bosonic NRG can likely be 
found in the rapidly developing field of ultracold bosonic 
gases, where indeed various realizations of spin-boson 
physics have been proposed 51 . 

Further, the physics of decohcrcncc of qubits naturally 
leads to variants of the spin-boson model. Interestingly, 
the description of 1 // noise in electrical circuits leads to 
sub-Ohmic damping with s — (at least over a certain 
range of energies). In this sub-Ohmic parameter regime, 
the bosonic NRG is one of the few methods which can 
give reliable answers, including, e.g., the existence of a 
quantum phase transition for < s < 1 - note that this 
transition does not appear in the popular non-interacting 
blip (NIB A) approximation 15 . Other modifications of 
standard spin-boson physics include the influence of lo- 
calized modes which interact with the qubit of interest - 
those modes can be represented by a discrete spin system, 
leading to so-called central spin models 17 . Usually, such 
systems map onto spin-boson models with a spectral den- 
sity consisting of a continuous (e.g. Ohmic) background 
and sharp peaks at certain frequencies 16 ; these models 
can be easily studied using NRG. 
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APPENDIX A: CALCULATION OF THE 
PARAMETERS OF THE SEMI-INFINITE CHAIN 

In this appendix, we describe the orthogonal trans- 
formation from the star-Hamiltonian (11) to the semi- 
infinite chain form [the chain-Hamiltonian (14)] and 
present equations relating the parameters of the two 
Hamiltonians. 

We start from the star-Hamiltonian (11): 

oo oo 

H s = H\ oc + 

n=0 Z V n n=0 

Our goal is to transform it to a semi-infinite chain 
eq. (14): 

H c = H loc + ] ^^ (bo + bl) 

OO 

+ [ e n b l h n + tn (&J,&n+l + b l+l b n)] • (A2) 

n=0 
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Here the main difference to the fermionic case (such as 
the Kondo model studied in Ref. f) is that the bosonic 
spectral function J(lu) is restricted to positive frequencies 
only. This asymmetry of J(co) influences the structure of 
the semi-infinite chain Hamiltonian (additional on-site 
energies e„ appear in eq. (A2) which are not present for 
particle- hole symmetry in the fermionic case). 

In the following we define a real orthogonal transfor- 
mation U, 



bn ^ ^ U nrn d m , 



(A3) 



m=0 



with U T U = UU T = 1, U* = U, so that the inverse 
transformation reads 



(A4) 



m=0 



Comparing the coupling terms between spin and bosons 
in H s and H c gives 



bo = —= } j n a n ■ 



so that 



U, 



On 



In 



(A5) 



(A6) 



The bosonic commutation relation [60, b ] = 1 applied to 
eq. (A5) gives 



Vo 



n=0 



J(x) dx . 



(A7) 



We are left with the equivalence of the free bosonic part 
in H s and H c : 



00 00 

£« a n a « = E [ £nb n bn + *« + b i +1 b r 



n=0 



n=0 



(A8) 



To obtain the recursion relations for e n and t n , we first 
put eq. (A4) into the left-hand side of eq. (A8) (for the 
annihilation operators only). We then sort the resulting 
equation for the operators b m . Comparing the prefactors 
of the terms containing b m we obtain for the operator bo'- 



X! £nOn U 0n = £0&o + *0&I , 



(A9) 



n=0 

and for b m with m > 0: 



trnb m+ i + t m —lb m _ 1 . (A10) 



n=0 



The expression for eo can be obtained from taking the 
commutator between bo and eq. (A9): 



e — ^ &i^0n i 
n=0 

From eq. (A9), we also obtain 

00 

tob\ = ^{£,n- eo) U 0n a n , 

n=0 

which gives immediately 

Uln — — (£,n - £o) Uo n . 

to 



(All) 



(A12) 



(A13) 



The value of to can be calculated by taking the commu- 
tator with the corresponding adjoint operator on both 
sides of eq. (A12). This results in 



to 



1 

rjo 



E«» 



\2 2 

eo) 7n 



Ln=0 



(A14) 



Equations (A6), (All), (A13), and (A14) initialize the 
recursion relations for the calculation of e m , t m and U mn . 
These recursion relations can be obtained by starting 
with eq. (A10) and proceeding in a similar way as above. 
The commutator beween b m and eq. (A10) gives 



= E^ 



2 

mn ' 



(A15) 



n=0 



From eq. (A10) we also find 



tmb m +i — ^ ^ {£,nUmn e m t> ^ ran t m — \Um— l n ) Cl^ . 



n=0 



(A16) 

From this equation, we obtain the expression for U m+ i n : 



U m +ln — " [(£n e m ) U mn t m — lU m — \ n ] . (A17) 
tm 

The values of t m can be calculated by taking the com- 
mutator with the corresponding adjoint operator on both 
sides of eq. (A16). This results in 



tm — 



^ y [i^n e m ) Umn tm— lU m — In] 



n=0 



(A18) 



Equations (A15), (A17), and (A18) complete the recur- 
sion relations for the calculation of the parameters of the 
chain Hamiltonian (14). 

Despite the simple structure of the input spectral func- 
tion, J(ui) = 2-Kaus 3 [s > 0), we did not succeed in solv- 
ing the recursion relations analytically (this is in fact 
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possible for the particle-hole symmetric soft-gap Ander- 
son model, where the hybridization function vanishes at 
the Fermi level as A(u>) = A |c<j| r , see Ref. 32; due to 
the particle-hole symmetry, the e„ vanish and the recur- 
sion relations have a much simpler structure) . Therefore 
the recursion relations have to be iterated numerically, 
in a similar way as for the fermionic case. Note that 
the derivation of the chain-Hamiltonian in the asymmet- 
ric fermionic case, where e n 7^ 0, is very similar to the 
bosonic case described above. The only differences are 
the structure of the coupling between impurity and the 
bath, and the fact that all commutators have to be re- 
placed by anticommutators. For a recent application of 
the NRG to a fermionic model with an asymmetric hy- 
bridization function, see Ref. 52. 

The resulting parameters of the chain-Hamiltonian, e n 
and t n , both fall off as A~"; in contrast to the fermionic 
case where t n oc A~"/ 2 . For large n, the ratio t n /e n 
approaches an s-dependent value. 



APPENDIX B: OPTIMAL BOSONIC BASIS IN 
THE STAR-NRG 

In this appendix, we present details of how we imple- 
ment the optimal basis for the bosons in the star-NRG 
to overcome the problem of the boson number divergence 
when the flow is to the localized fixed point (see the dis- 
cussions in Sees. HE and HID). 

In each step of the star-NRG, a new bosonic degree of 
freedom is added to the Hamiltonian. The rcnormaliza- 
tion group transformation is given by eq. (19): 



Hn+1,8 

+ A N+1 



A-ff/v,s 

Civ+iajv+i a JV+i + 2^^ N+1 ( aN + 1 + a lv+i) 

(Bl) 



The problem discussed in Sec. HE is that upon ap- 
proaching the localized fixed point in the sub-Ohmic case, 
the displacements 9 m for the bosonic site N increase ex- 
ponentially with N, see eq. (36). The displacements can 
be taken into account by constructing an appropriate ba- 
sis \s(N + 1)) for the new bosonic degree of freedom. 

To construct this basis, we start with a simplified 
Hamiltonian of the form: 



H = a)a + 9o z (a f + a) 



(B2) 



and proceed as follows: in the first step, we set up an 
optimized basis for H (optimized in the sense that the 
lowest lying eigenstates of H are described with only a 
small number of basis states). Then we use this basis, 
denoted as \s(N + l))g, for the actual NRG iteration, 
and finally, we describe a self-consistent procedure to de- 
termine the parameter 9. 

Let us first discuss how to construct the optimized ba- 
sis for H. Consider the following operators: 



H± g =a*a±6 (a 1 * + a) 



(B3) 



The eigenstates of H±g are denoted as \m)±$ (m = 
0, 1, ...). We obtain 



H± e \m)± e = (to - 9 2 ) \m)± e ■ 



and 



\m)±e 



( at - a )|m) 



(B4) 



(B5) 



with I to) the eigenstates of a^a. The basis states should 
describe the +9 and —9 displacements on an equal foot- 
ing; therefore we proceed with symmetrized eigenstates 
l m ) c /o constructed in the following way: 



He = c c . m [\m)g + (-l) m |TO>-e] 
|to) = c Q . m [\m)g - (-l) m |TO,)_ e ] 

n 1 N * 1 

TO = 0, 1, — - 1 , 



(B6) 



with normalization constants c e / 0iTO . Note that here we 
have to choose an even number Ny,. The even and odd 
parity states are orthogonal to each other, e (n|TO) = 0, 
whereas states with the same parity are not necessar- 
ily orthogonal. An orthogonalization procedure for both 
even and odd parity states then gives the final set of basis 
states: 



|0)e 
|l)c 
|2)e 



|0)c 

Ce,l{|l>e 

Cc, 2 {|2) c 



3<0|l)e|0)e} 
a(l|2) e |l)e- 



3<0|2>e|6> e } 



(B7) 



with normalization constants C c / o m . The same orthogo- 
nalization is performed for the odd parity states. In this 
way, we obtain orthogonal states, characterized by 
the parameter 9, which form the basis \s(N + l))$ for the 
diagonalization of £fjv+i, s : 

\ S (N + l))g = {|0} e , |l) c , . . . , |0) o , |1) , . . .} . (B8) 

The calculation of the matrix elements HN+i, s {rs,r' s') 
(see eqs. (43-46)) involves matrix elements of the form 



e (s{N + l)\ a N+1 + a N+1 \s'(N + 



(B9) 



To evaluate these matrix elements [and the scalar prod- 
ucts in cq. (B7)] we have to express the states \fh) e / 

in terms of the eigenstates \n) of o] v+1 ajv+i. This can 
be performed using the following recursion relations for 

(n\m) g : 



(n\m + l)e 
(n\0)g 
(0\m) e 



9 



Vto + 1 



(n\m) e + 



VmTT 



(- 



'to! 



(n- l\m) e , 



(BIO) 
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(and for (n|m)_g by replacing 9 by —9). The summation 
over n in the calculation of matrix elements and scalar 
products has to be performed numerically which limits 
the number of states \n) to some finite, although very 
large, value L (values up to L w 10 7 can be used). To 
construct an optimized basis for the displacements 9, L 
should be large enough (at least of the order of 9 2 ) to 
include a sufficient number of states |n) in the calculation. 

For the special case of A = 0, the parameter 9 for the 
construction of the basis \s(N+l))g is exactly known (see 
Sec. HE). This is different in the general case of finite A 
where we have to find a scheme to determine the optimal 
value 9* . The general strategy to find this optimal value 
has been discussed in Sec. II F. For the actual numerical 
calculation it turns out that the following self-consistent 
scheme is much more efficient. 

For the ground state \g) of H (B2) the expectation 
value (g\a^a\g) is equal to 9 2 . We use this relation to 
determine the 9 used for the NRG calculation: 



JV+l 



(s|aj\r+i a iv+i|ff)Ar+i 



(Bll) 



where Iff) jv+i is the ground state of _ff/v + i. s which 
has been obtained from diagonalizing the matrix 
HN+i,s(rs, r's') using the basis \s(N + l))$ characterized 
by the parameter 9. In other words, eq. (Bll) defines a 
self-consistent scheme to calculate 9 for each NRG step. 

The converged value 9* gives the optimal basis for 
adding the site N + 1 in the NRG iteration. It corre- 
sponds to the value 9* which characterizes the minimum 
of the energy levels in Fig. 4. The energy levels calcu- 
lated in this way show a much weaker dependence on N\> 
which leads, for example, to the rapid convergence of rib 
with increasing N\, as shown in Fig. 13. 

After the diagonalization of i?jv+i, s with the optimized 
basis \s(N+l))g* , we can continue the NRG iteration by 
adding the site N + 2. 
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